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ABSTRACT 

We introduce a simple model to self-consistently connect the growth of galaxies to 
the formation history of their host dark matter halos. Our model is defined by two 
simple functions: the "baryonic growth function" which controls the rate at which 
new baryonic material is made available for star formation, and the "physics function" 
which controls the efficiency with which this material is converted into stars. Using 
simple, phenomenologically motivated forms for both functions that depend only on 
a single halo property, we demonstrate the model's ability to reproduce the z=0 red 
and blue stellar mass functions. Furthermore, by adding redshift as a second input 
variable to the physics function we show that the reproduction of the global stellar 
mass function out to z=3 is improved. We conclude by discussing the general utility 
of our new model, highlighting its usefulness for creating mock galaxy samples which 
have a number of key advantages over those generated by other techniques. 

Key words: galaxies: evolution - galaxies: formation - galaxies: halos - galaxies: 
statistics - galaxies: stellar content 
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1 INTRODUCTION 

Theoretical models are an important and commonly used 
tool for interpreting and furthering our understanding of ob- 
served galaxy populations. Typically, these models are used 
to generate mock galaxy catalogues that can be compared to 
equivalent samples drawn from the real Universe. Knowledge 
of the models' construction, combined with their successes 
and failures in reproducing the observations, can often allow 
important inferences to be made about the physics of galaxy 
formation and evolution. 

There are a number of different methods for generat- 
ing mock galaxy samples for comparison with observations. 
At the most advanced and complex end of the scale are 
full hydrodynamic simulations. These attempt to solve the 
physics of galaxy formation from first principles, directly 
modelling complex baryonic processes such as cooling and 
shocks in tandem with the dissipationless growth of dark 
matter structure. Unfortunately, the associated high compu- 
tational cost prohibits the resolution of small scale physical 
processes such as star formation and black hole feedback in a 
volume large enough to provide a cosmologically significant 
sample of galaxies. Hence hydrodynamic simulations often 
resort to parametrised approximations to deal with these 
unresolved "sub-grid" processes. In addition they must also 
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deal with the complex and often poorly understood numer- 
ical effects that come with modelling dissipational physics 
using finite physical and temporal resolutions. 

Semi-analytic galaxy formation models attempt to over- 
come the computational costs associated with hydrodynamic 
simulations by separating the baryonic physics of galaxy for- 
mation from the dark matter dominated growth of struc- 
ture (White & Frenk 1991; Kauffmann et al. 1999). This is 
achieved by taking pre-generated dark matter halo merger 
trees and post-processing them with a series of physically 
motivated parametrisations that attempt to capture the 
mean behaviour of the dominant baryonic processes involved 
in galaxy formation. The resulting speed means that these 
models can be used to generate cosmologically significant 
samples of galaxies using only modest computing resources. 
However, semi-analytic models typically require a number 
of free parameters, many of which are often not well con- 
strained by theory or observation (Neistein & Weinmann 
2010). The complicated and degenerate nature of the dif- 
ferent physical prescriptions also means that the effects of 
these parameters on the final galaxy population are often 
highly degenerate and can be difficult to interpret (e.g. Lu 
et al. 2012). Additionally, our relatively poor understanding 
of high redshift galaxy formation means that at least some of 
the parametrisations used may not be appropriate at these 
early times (Mutch et al. 2013; Henriques et al. 2012). 

For many science questions and applications we are not 
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required (or able) to include and understand all of the rel- 
evant input physics. In these cases it is often sufficient to 
construct simple "toy" models (e.g. Tacchella et al. 2012; 
Dekel et al. 2013). These typically build an average popula- 
tion of galaxies using vastly simplified approximations, and 
are designed to test new ideas and interpretations or to al- 
low the investigation of particular trends or features found 
in observational data. One example is the "reservoir" model 
of Bouche et al. (2010). Here, averaged dark matter halo 
growth histories are used to track the typical build up of 
cold gas in galaxies. In their fiducial formalism, the accre- 
tion of baryons onto a galaxy is only allowed to occur when 
the host halo lies in a fixed mass range. However, within 
this range the accretion is modelled as a simple fraction of 
the halo growth rate. Using a standard Kennicutt-Schmidt 
law (Kennicutt 1998) for star formation, this simple model 
is able to reproduce the observed scaling behaviours of the 
star forming main sequence and Tully-Fisher relations. Ex- 
panding upon this framework, Krumholz & Dekel (2012) 
also introduced a metallicity dependent star formation ef- 
ficiency, allowing them to straightforwardly investigate the 
associated effects on the star formation histories of galaxies. 

Rather than attempting to generate galaxy populations 
based on our theories of the relevant physics, an alternative 
method of generating mock galaxy samples is to use purely 
statistical methods. Halo Occupation Distribution (HOD) 
models use observed galaxy clustering measurements to con- 
strain the number of galaxies of a particular type within a 
dark matter halo of a given mass (Peacock & Smith 2000; 
Zheng et al. 2005). For the purpose of constructing mock 
galaxy catalogues, such a methodology has the advantage 
that it requires no knowledge of the how each galaxy forms 
and is also statistically constrained to produce the correct 
result. However, this limits our ability to learn about the 
physics of galaxy evolution, as one has no way to self- 
consistently connect individual galaxies at any given redshift 
to their progenitors or descendants at other times. 

A similar method to HODs for creating purely sta- 
tistical mock catalogues is sub-halo abundance matching 
(SHAM; e.g. Conroy et al. 2006). SHAM models are typi- 
cally constructed by generating a sample of galaxies of vary- 
ing masses, drawn from an observationally determined stel- 
lar mass function. Each galaxy is then assigned to a dark 
matter halo taken from a halo mass function generated us- 
ing an N-body simulation. This assignment is made such 
that the most massive galaxy is placed in the most massive 
halo and so on, proceeding to lower and lower mass galax- 
ies. Due to possible differences in the formation histories of 
halos of any given mass, an artificial scatter is often added 
during this assignment procedure (e.g. Conroy et al. 2007; 
Behroozi et al. 2012; Moster et al. 2013). 

By leveraging the use of dark matter merger trees as 
the source of the halo samples at each redshift, both HOD 
and SHAM studies have been able to provide important con- 
straints on the average build up of stellar mass in the Uni- 
verse (Zheng et al. 2007; Conroy & Wechsler 2009; Moster 
et al. 2013; Behroozi et al. 2012). This has allowed these 
studies to also draw valuable conclusions about processes 
such as the efficiency of star formation as a function of halo 
mass and the role of intra-cluster light (ICL) in our account- 
ing of the stellar mass content of galaxies. However, both 
HOD and SHAM models are applied independently at indi- 



vidual redshifts and do not self-consistently track the growth 
history of individual galaxies. This limits the remit of these 
models to considering only the averaged evolution of certain 
properties over large samples. 

Our goal in this work is to present an alternative class of 
galaxy formation model which allows us to achieve the "best 
of both worlds" ; providing both a self-consistent growth his- 
tory of each individual galaxy, whilst also minimising any as- 
sumptions about the physics which drives this growth. This 
is achieved by tying star formation (and hence the growth 
of stellar mass) to the growth of the host dark matter halo 
in N-body dark matter merger trees using a simple but well 
motivated parametrisation that depends only on the prop- 
erties of the halo itself. In this way, we are able to provide a 
complete formation history for every galaxy. The model we 
present is closely related to that of Cattaneo et al. (2011), 
but with a number of important generalisations that increase 
its utility whilst still maintaining a high level of transparency 
and simplicity. 

This paper is laid out as follows: In §2 we introduce the 
framework of our new model. In particular, §2.2 focusses on 
how we build up the baryonic content of dark matter halos, 
with the practical details of the model's application outlined 
in §2.3. In §3 we present some basic results, in particular 
investigating the model's ability to reproduce the observed 
galaxy stellar mass function at multiple redshifts. In §4 we 
discuss our findings as well as outline the general utility of 
the model and a number of possible ways in which it can be 
extended. Finally, we present a summary of our conclusions 
in §5. 

A WMAPl (Spergel et al. 2003) ACDM cosmology with 
r2m=0.25, Qa=0.75, ^^6=0.045 is utilised throughout this 
work. In order to ease comparison with the observational 
datasets employed, all results are quoted with a Hubble con- 
stant of h—^.7 (where /i=ifo/100kms~"^Mpc~"^) unless oth- 
erwise indicated. Magnitudes are presented using the Vega 
photometric system and a standard Salpeter (1955) IMF is 
assumed throughout. 



2 THE SIMPLEST MODEL OF GALAXY 
FORMATION 

2.1 The growth of structure 

The aim of the model presented in this work is to self- 
consistently tie the growth of galaxy stellar mass to that 
of the host dark matter halos in as simple a way as possible. 
In order to achieve this we require knowledge of the proper- 
ties and associated histories of a large sample of dark matter 
halos spanning the full breadth of cosmic history. We obtain 
this in the form of merger trees constructed from the output 
of the N-body dark matter Millennium Simulation (Springel 
et al. 2005). 

Using the evolution of over 10^*^ particles in a cubic 
volume with a side length of 714 Mpc, the Millennium Sim- 
ulation merger trees track the build up of dark matter halos 
larger than approximately 2.9xlO^'^M0 over 64 temporal 
snapshots. These snapshots are logarithmically spaced in ex- 
pansion factor between redshifts 127 and 0, with an average 
separation of around 200 Myr. Each individual dark mat- 
ter structure is identified using a Friends-of- Friends linking 
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algorithm with further substructures (sub-halos) identified 
using the SUBFIND algorithm of (Springel et al. 2001). The 
simulation employs a concordance ACDM cosmology com- 
patible with first year WMAP (Spergel et al. 2003) param- 
eters: (Qm,^A,cr8,/io) = (0.25,0.75,0.9,0.73). 

2.2 The baryonic content of dark matter halos 

The maximum star formation rate of a galaxy is regulated 
by the availability of baryonic material that can act as fuel. 
In our formation history model we assume that every dark 
matter halo carries with it the universal fraction of baryonic 
material, /b=0.17 (Spergel et al. 2003). However, some of 
these baryons will already be locked up in stars or contained 
in reservoirs of material that are unable to participate in 
star formation. Therefore we parametrise the dependence 
of the amount of newly accreted baryonic material which 
is available for star formation on the properties of the host 
dark matter halo using a baryonic growth function, Fgrowth • 

In practice, only some fraction of this available mate- 
rial will actually make its way in to the galaxy, with an 
even smaller amount then successfully condensing to form 
stars in a suitably short time interval. The efficiency with 
which this occurs depends on a complex interplay of non- 
conservative baryonic processes, both internal and external 
to the galaxy-halo system. A number of important exam- 
ples include shock heating, feedback from supernova and 
active galactic nuclei (AGN), as well as environmental pro- 
cesses such as galaxy mergers and tidal stripping. Here we 
assume that all of these complicated and intertwined mecha- 
nisms can be distilled down into a single, arbitrarily complex 
physics function, Fphys- 

Combining all of this together, we can write down a 
deceptively simple equation for the growth rate of stellar 
mass (M*) in the universe on a per-halo basis: 



Lookback time (Gyr) 
7.9 10.6 11.9 12.5 12.913.2 



M* 



- phys ■ 



(1) 



In the following sections, we discuss the form we employ 
for the baryonic growth and physics functions in turn. 

2.2.1 The baryonic growth function 

In order to explore the simplest form of our formation his- 
tory model, we begin by assuming that as a dark matter halo 
grows, all of the fresh baryonic material it brings with it is 
immediately available for star formation. This corresponds 
to a baryonic growth function which is simply given by the 
rate of growth of the host dark matter halo: 



growth 



/b 



dMvir 

dt 



(2) 



In practice halos of the same z=0 mass may show a 
diverse range of growth histories, all of which are captured 
by our model. In Fig. 1 we demonstrate this by showing 
the individual growth histories of a random sample of dark 
matter halos selected from the Millennium simulation in two 
narrow mass bins. From this figure we see that there can be 
significant variations in the time at which similar halos at 
redshift zero reach a given mass. For example, in the up- 
per halo mass sample, some halos reach 10^^ M© by z=5 
whilst others do not reach this value until z=2. In addition, 
some halos may have complex growth histories, achieving 




2 3 

Redshift 

Figure 1. This figure demonstrates the large variation in the 
possible growth histories of halos which all have approximately 
equal masses by redshift zero. The blue and red lines represent 30 
randomly selected growth histories for halos with final Mvir val- 
ues of approximately IQ-*^^ and IQ-'^'^ Mq respectively. Variations 
of 3-4 Gyr in the time at which these halos reach a given mass 
is common. Unlike statistical techniques for tying galaxy prop- 
erties to their host halos, our formation history model implicitly 
includes the full range of different halo growth histories and their 
effects on the predicted galaxy population. 



their maximum mass at z>0. This can potentially be caused 
by a number of processes such as stripping during dynami- 
cal encounters with other halos. Since the baryonic growth 
function maps the formation history of each individual dark 
matter halo to the stellar mass growth of its galaxy, this 
diversity in growth histories is fully captured, propagating 
through to be refiected in the predicted galaxy populations 
at all redshifts. This is an important attribute of our model 
that sets it apart from other statistical-based methods which 
merely map the properties of galaxies to the instantaneous 
or mean properties of halos, independently of their histories 
(e.g. HOD and SHAM models). These methods typically 
have to add artificial scatter to approximate the effects of 
variations in the halo histories, whereas this variation is a 
self-consistent input to our formation history model. 



2.2.2 The physics function 

The physics function describes the efficiency with which 
baryons are converted into stars in halos of a given mass. 
The form of this function may be arbitrarily complex, how- 
ever the goal of this work is to find the simplest model which 
successfully ties the growth of galaxy stellar mass to the 
properties of the host dark matter halos. The physics func- 
tion is not meant to provide an accurate reproduction of the 
details of the full input physics, but rather their combined 
effects on the growth of stellar mass in the Universe. In this 
spirit, we begin by assuming that there is only one input 
variable: the instantaneous virial mass of the halo, Mvir- 

Although still not understood in detail, the observed 
relationship between dark matter halo mass and galaxy stel- 
lar mass is well documented (e.g. Zheng et al. 2007; Wang 
et al. 2012; Yang et al. 2012). Assuming the favoured ACDM 
cosmology, a comparison of the observationally determined 
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galactic stellar mass function to the theoretically determined 
halo mass function indicates that the averaged efficiency of 
stellar mass growth varies strongly as a function of halo 
mass. In Fig. 2, we contrast a Schechter function fit of the 
observed redshift zero stellar mass function (Bell et al. 2003, 
solid blue line) against the dark matter halo mass function of 
the Millennium Simulation (red dashed line). The halo mass 
function has been multiplied by /b in order to approximate 
the total amount of baryons available for star formation in 
a halo of any given mass. 

The increased discrepancy between the stellar mass 
function and halo mass functions at both low and high 
masses indicates that the efficiency of star formation is 
reduced in these regimes. It is commonly held that at 
low masses the shallow gravitational potential provided by 
the dark matter halos allows supernova feedback to effi- 
ciently eject gas and dust from the galaxy. This reduces 
the availability of this material to fuel further star forma- 
tion episodes, hence temporarily stalling in-situ stellar mass 
growth. Other processes such as the photo-ionisation heat- 
ing of the intergalactic medium may also play an important 
role in reducing the efficiency of star formation in this low 
mass regime (Benson et al. 2002, and references therein). 
At high halo masses, it is thought that inefficient cooling 
coupled with strong central black hole feedback also leads 
to a quenching of star formation (e.g. Croton et al. 2006). 
Therefore, it is only between these two extremes, around the 
knee of the galactic stellar mass function, that stellar mass 
growth reaches its highest average efficiency. 

We begin by parametrising the physics function as a 
simple log-normal distribution centred around a halo virial 
mass Mpeak, and with a standard deviation cfm^-^^'- 



i^phys(Mvir/M0) = Em^;, exp 



f AM^ 



\ O-M^i^ 



(3) 



where AMvir=logio(Mvir/M0)- log^o (Mpeak /M©) and the 
parameter Sm^^^ represents the maximum possible efficiency 
for converting in- falling baryonic material into stellar mass, 
achieved when Mvir=Mpeak- Such a distribution was found 
by the SHAM study of Conroy & Wechsler (2009) to pro- 
vide a good match to the derived star formation rates as a 
function of halo mass for z<l. 

This simple form of the physics function provides a 
number of desirable properties. In Fig. 3, we present the 
average growth histories of five samples of dark matter ha- 
los chosen from the Millennium Simulation merger trees by 
their final redshift zero masses (solid blue lines). For clar- 
ity, we only plot these histories out to redshifts where more 
than 80% of the halos in each sample have masses which 
are twice the resolution limit of the input merger trees. The 
grey shaded region indicates the amplitude of the physics 
function defined by Eqn 3 when using our fiducial param- 
eter values (see §3.1 for details). As the halos grow, they 
pass through the region of efficient star formation at differ- 
ent times depending on their final masses. Galaxies hosted 
by the most massive z—0 halos form the majority of their in- 
situ stellar mass at earlier times whereas those in the lowest 
mass halos are still to reach the peak of their growth. In ad- 
dition, lower mass halos tend to spend a longer time in the 
efficient star forming regime compared to their high mass 
counterparts. These trends qualitatively agree with the ob- 
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Figure 2. A comparison of the observed galactic stellar mass 
function (blue solid line) and the halo mass function of the Mil- 
lennium Simulation (red dashed line). The halo mass function 
has been multiplied by the universal baryon fraction in order to 
demonstrate the maximum possible stellar mass content as a func- 
tion of halo mass. The closer the stellar mass function is to this 
line, the more efficient star formation is in halos of the correspond- 
ing mass. If galaxies were to form stars with a fixed efficiency at 
all halo masses then the slope of the stellar mass function would 
be identical to that of the halo mass function. The differing slopes 
at both high and low masses indicates that star formation (as a 
function of halo mass) is less efficient in these regimes. At low 
masses, this is commonly attributed to efficient gas ejection due 
to supernova feedback, whereas at high masses energy injection 
from central super-massive black holes is thought to be able to 
effectively reduce the efficiency of gas cooling. However, many 
other physical processes may also contribute in both regimes. 



served phenomenon of galaxy downsizing (e.g. Cowie et al. 
1996; Cattaneo et al. 2008). 

Sub-halo abundance matching studies have suggested 
that Vmax may be more tightly coupled to the stellar mass 
growth of galaxies than Mvir (e.g. Reddick et al. 2012). This 
makes intuitive sense as Vmax is directly related to the gravi- 
tational potential of the inner regions of the host halo, where 
galaxy formation occurs. Therefore, as well as virial mass 
we also consider the case of a physics function where the 
dependent variable is the instantaneous maximum circular 
velocity of the host halo, Vmax: 



i^phys(Kiax/(kmS ^)) = ^Vmax ^Xp 






(4) 



where AKiax=logio(Knax/(kms~^))-logio(V^peak/(kms~^)). 
To avoid confusion, from now on we will refer to the forma- 
tion history model constructed using this physics function 
as the "static Vmax model". Similarly, we will refer to the 
case of Fphys(Mvir) as the "static Mvir model". 

In Fig. 3 we show the average Vmax growth histories 
for a number of different z—0 selected samples. The y-axis 
has been scaled such that the grey band also correctly de- 
picts the changing amplitude of the Vmax physics function as 
well as its Mvir counterpart. Additionally, each of the Vmax 
samples in Fig. 3 (red dashed lines) are chosen to have mean 
z—0 values close to that of the five Mvir samples (blue lines). 
However, there are clear differences between the growth his- 
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Figure 3. The mean virial mass (Mvir) growth histories for five samples of dark matter halos with varying final masses (blue solid lines). 
Each sample is only plotted out to a redshift limit determined by where 80% of the halos contain more than 40 particles. The grey shaded 
region indicates the amplitude of the physics function in Eqn. 3. This is also illustrated by the log-normal curve on the right-hand side 
(orange shaded region). Galaxies with halo masses within the peaked region form stars efficiently; outside this mass range, the amount 
of star formation is negligible. All galaxies of sufficient mass at z=0 will cross the efficient star formation mass range at some point in 
their history, with this period typically coming earlier for more massive halos. Also shown is the mean maximum circular velocity (T/max) 
evolution for halos of varying final velocity (red dashed lines). These samples have been chosen to have z=0 maximum circular velocity 
values similar to the mean values of the five mass-selected samples. There are clear differences in the evolution of Vmax and Mvir which 
will in turn result in differences between the produced galaxy populations. 



tories of these two halo properties. In particular, the evolu- 
tion of Vmax is slightly flatter, resulting in halos transition- 
ing out of the efficient star forming region at an earlier time 
than the equivalent Mvir sample. Such differences will have 
important consequences for the time evolution of the galaxy 
populations generated by each of the two physics functions 
and we highlight some of these in §3 below. 

By combining the baryonic growth function with a 
physics function of the forms presented here, our resulting 
model may be thought of as a simplified and extended ver- 
sion of that presented by Bouche et al. (2010). Unlike their 
model, the scaling of gas accretion efficiency with halo mass, 
and the dependence of star formation on previously accreted 
material, is implicitly contained within our physics function. 
Most importantly though, Bouche et al. (2010) use statisti- 
cally generated halo growth histories instead of simulated 
merger trees. Hence their model contains no information 
about the scatter due to variations in halo formation his- 
tories as well as no self- consistent stellar mass growth due 
to mergers. 

The model of Cattaneo et al. (2011) also uses simulated 
merger trees as input and thus shares many of the same 
advantages as our formation history model. However, their 
model ties the properties of galaxies to the instantaneous 



properties of their host halos alone. In contrast, we use the 
full information of the mass accretion history to describe the 
availability of baryonic material for star formation. Also, we 
make no attempt to motivate the precise form of our model 
in terms of combinations of particular physical processes and 
their scalings with halo properties, as Cattaneo et al. (2011) 
do. This allows our model to remain maximally general and 
flexible. 



2.3 Generating the galaxy population 

Armed with the forms of our baryonic growth function 
(Eqn. 2) and physics function (Eqn. 3), we now discuss the 
practical implementation of the formation history model to 
generate a galaxy population from the input dark matter 
merger trees. 

For each halo in the tree, the change in dark matter halo 
mass, coupled with the time between each merger tree snap- 
shot, provides us with the value of dM^ir/dt. This change in 
mass, naturally includes growth due to both smooth accre- 
tion and merger events. Combined with the instantaneous 
value of Mvir or Vmax we can calculate a star formation rate 
for the occupying galaxy following Eqn. 1. 

Some fraction of the mass formed by each new star for- 
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mation episode will be contained within massive stars. The 
lives of these stars will be relatively short and therefore they 
will not contribute to the measured total stellar mass con- 
tent of the galaxy. In order to model this effect we invoke the 
"instantaneous recycling" approximation (Cole et al. 2000), 
whereby some fraction of the mass of newly formed stars is 
assumed to be instantly returned to the galaxy ISM. Based 
on a Salpeter (1955) initial mass function we take this frac- 
tion to be 30%, however, we note that changes to this value 
can be trivially taken into account by appropriately scaling 
the value of 8 in the physics function. 

Although well motivated and conceptually simple, our 
use of dM^ir/dt in the baryonic growth function (Eqn. 2) 
does introduce some practical considerations. For example, 
the change in halo mass from snapshot-to-snapshot in the 
input dark matter merger trees can be stochastic in nature, 
especially for the case of low mass or diffuse halos identified 
in regions of high density. Also, when satellite galaxies fall 
into larger systems their halos are tidally stripped, leading 
to a negative change in halo mass and thus a reduction in 
stellar mass according to Eqn. 1. In the real universe, we 
expect that the galaxy is located deep within the poten- 
tial well of its host halo and is therefore largely protected 
from the earliest stripping effects suffered by the dark mat- 
ter (Pefiarrubia et al. 2010). We must therefore decide when, 
if at all, to allow stellar mass loss when using this formalism. 
For simplicity, we address this by setting the star formation 
rate of satellite galaxies to be zero at all times; in other 
words fixing their stellar mass upon in fall. This is unlikely 
to be true in the real universe across all mass and environ- 
ment scales (Weinmann et al. 2006), however, the assump- 
tion of little or no star formation in satellite galaxies is a 
reasonable approximation and is relatively common in ana- 
lytic galaxy formation models (e.g. Kauffmann et al. 1999; 
Cole et al. 2000; Bower et al. 2006; Croton et al. 2006). It is 
also in keeping with our goal of finding the simplest possible 
model. 

The form of the baryonic growth function presented in 
Eqn. 2 above is only one of a number of possibilities. As an 
example, one could use the instantaneous halo mass divided 
by its dynamical time, Mvir/td^n- This quantity grows more 
smoothly over the lifetime of a halo and is never negative. 
Additionally, one may speculate that this is a better rep- 
resentation of the link between stellar and halo mass build 
up. However, for simplicity, we do not investigate alterna- 
tive forms of the baryonic growth function, but leave this to 
future work. 

Satellite galaxies are explicitly tracked in the input 
merger trees until their host sub-halos can no longer be iden- 
tified or fall below the imposed resolution limit of 20 par- 
ticles. At this point, their position is approximated by the 
location of the most bound particle at the last snapshot the 
halo was identified. We then follow Croton et al. (2006) in 
assuming that the associated satellite galaxy will merge with 
the central galaxy of the parent halo after a timescale moti- 
vated by dynamical friction arguments (Binney & Tremaine 
2008): 



tn 



1.17 



Kirr,' 



G rrisat ln(l + Mvir/msat) 



(5) 



where V^ir and Mvir are the virial velocity and mass of 
the parent dark matter halo in kms~^ and M© respec- 



tively, Tsat is the current radius of the satellite halo in 
kpc, and msat is the mass of the satellite in Mq. In these 
units, the gravitational constant, G, is given by 4.40 x 
10"^ kpc^ km s-^ Mq^ MyY-\ 

The ensuing merger event results in a galaxy with a stel- 
lar mass equal to the direct sum of the progenitor masses. 
Also, as in- falling satellite halos are stripped of their dark 
matter this mass is typically added to the parent halo, con- 
tributing to its dMvir/dt value and thus the amount of star 
formation in the central galaxy. Combined with our simple 
baryonic growth function that assumes all of the incoming 
baryonic material is available for star formation (irrespec- 
tive of whether or not it is already locked up in stars), our 
model implicitly includes merger driven star-bursts with an 
increased efficiency. 

Knowledge of the star formation rate of each galaxy also 
allows us to calculate luminosities. For this purpose we use 
the simple stellar population models of Bruzual & Chariot 
(2003) and assume a Salpeter (1955) initial mass function. In 
the real universe supernova eject a enriches the intra galactic 
medium, altering the chemical composition of the next gen- 
eration of stars and the spectrum of the light they emit. As 
we do not track the amount of gas or metals in our simple 
model, we assume all stars are of 1/3 solar metallicity. This 
is a common assumption when no metallicity information is 
available. Finally, a simple "plane-parallel slab" dust model 
(Kauffmann et al. 1999) is applied to the luminosity of each 
galaxy in order to provide approximate dust extincted mag- 
nitudes. These magnitudes are used below to augment our 
analysis by allowing us to calculate the B—V colour for each 
galaxy at z=0. However, our main focus will remain on stel- 
lar masses as these are a direct model prediction. 



3 RESULTS 

Having outlined the methodology and implementation of our 
simple formation history model, we now present some basic 
results which showcase its ability to recreate observed distri- 
butions of galaxy properties. We begin by considering red- 
shift zero alone, before moving on to investigate the results 
at higher redshifts. Throughout, we contrast the variations 
between the predicted galaxy populations when using Mvir 
or Vmax as the dependant variable of the physics function 
(Eqns. 3 & 4). 



3.1 Redshift zero 

In order to determine the "best" parameter values for the 
Mvir model, we calibrate them against Schechter function 
fits of the observed red and blue galaxy stellar mass func- 
tions of Bell et al. (2003). This calibration was done using 
Markov Chain Monte-Carlo (MCMC) parameter estimation 
techniques (for details of our implementation see Mutch 
et al. 2013). The observed mass functions are constructed 
from a ^-band limited sample taken from a combination of 
Sloan Digital Sky Survey (SDSS) early release (Stoughton 
et al. 2002) and Two Micron All Sky Survey (2MASS; Jar- 
rett et al. 2000) data, with a magnitude dependent colour 
cut used to divide the red and blue galaxy populations. To 
similarly split the model galaxies into red and blue samples 
we employ a more basic mass/magnitude independent colour 
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Table 1. The fiducial parameter values of the physics function when using either Mvir or Vmax as the dependant variable. Values are 
presented for both the non-evolving (see §3.1) and evolving (see §3.3) form. Also shown are the ranges for the flat priors used during 
the MCMC calibration. The parameters Iog2o(^peak/M0) and log]^Q(ypeak/(kms~-'^)) indicate the halos which possess the peak star 
formation efliciency (at z=0 in the evolving case). The a and S parameters represent the width and height of the Gaussian physics 
function respectively. For the evolving models, a, /3 and 7 indicate the rate of power-law evolution of the peak location, width and height 
of the physics function respectively. The non-evolving model parameters were chosen to provide the best reproduction of the observed 
z=0 colour-split stellar mass function of Bell et al. (2003). In the evolving case, the values were chosen to additionally reproduce the 
evolution of the peak stellar-halo mass relation of Moster et al. (2013). 
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cut of B— y = 0.8. This is equivalent to the colour division 
found by the 2dF Galaxy Redshift Survey (Cole et al. 2005). 

For all of the results presented in this work we use a 
minimum of 130 000 model calls in the integration phase of 
our MCMC chains. Due to computational limitations we are 
unable to utilise the full Millennium simulation volume and 
instead restrict ourselves to a random sampling of 1/128 
of the total simulation merger trees. This is equivalent to 
a co-moving volume of approximately 9.8 x lO^Mpc^h"^. 
We use the same random merger tree set throughout. Flat 
priors were used for each parameter with ranges as presented 
in Table 1. To ensure that all of our MCMC chains are fully 
converged we employ the Rubin-Gelman statistic (Gelman 
Sz Rubin 1992) as well as visually inspect the chain traces. 




10.0 10.5 11.0 

logio(M,[M0]) 



12.0 



Using the posterior probability distributions of the 
MCMC fitting procedure allows us place 68 and 95% confi- 
dence limits on all of our model results, both constrained and 
predicted. In this work, confidence intervals are calculated 
from a large sample of model runs (60-200) that have pa- 
rameter combinations randomly sampled from the relevant 
posterior distributions. 

An important consideration when statistically calibrat- 
ing any model against observational data is the use of re- 
alistic observational uncertainties (Mutch et al. 2013). As 
discussed by Bell et al. (2003), there is likely significant 
systematic uncertainties associated with their stellar mass 
function estimation which are not formally included in the 
relevant Schechter function parameter values. To overcome 
this we utilise the uncertainties of Baldry et al. (2008) which 
are calculated by comparing the global mass functions that 
result from five independent stellar mass determinations of 
a single galaxy sample. We then partition this global un- 
certainty between red and blue galaxies. This is done such 
that the fractional contribution to the uncertainty due to red 
(blue) galaxies is equal to the fraction of red (blue) galaxies 
in each stellar mass bin. 



Figure 4. The colour-stellar mass diagram of the static Mvir 
model. The black dashed hue at B-V=^.d> (Cole et al. 2005) indi- 
cates the value used to divide the model galaxies into red and blue 
samples. There is a clear ridge of over-density extending across 
all stellar masses at B— y?^0.87 representing the red sequence. 



3.1.1 The Mvir model 

The best fit parameters for the static Mvir model are pre- 
sented in Table 1 (see §3.2 for the 'evolving' model). The 
preferred Mpeak value of 10^^"^ M© implies that galaxies 
in halos slightly less massive than that of the Milky Way 
(Mvir ~ lO"*^^ M0; Xue et al. 2008) are on average the most 
efficient star formers. In these halos, 56% of all freshly ac- 
creted baryonic material is converted into stars, as indicated 
by the value of Sm^^^ • 

In Fig. 4 we show the colour-stellar mass diagram pro- 
duced using the best parameters of our static Mvir model. 
The black dashed line indicates the colour split used to di- 
vide the galaxies into red and blue populations. Although 
there is a lack of a clear colour bi-modality as seen in obser- 
vational data (e.g. Baldry et al. 2004), we still find a clear 
overabundance of galaxies with B—V^O.'^l corresponding 
to the observed "red sequence" . The presence of this feature 
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Figure 5. The red and blue z=0 stellar mass functions produced by the static formation history model (solid lines) when using Mvir 
(left panel) and Vmax (right panel) as the input variable to the physics function (Eqns. 3 & 4). Galaxy colour is classified using a 
mass-independent colour cut of B—V=^.d>. The free model parameters have been calibrated to provide the best possible reproduction of 
the observed stellar mass function of Bell et al. (2003) (error-bars) and are presented in Table 1. Shaded regions indicate the 68% (dark) 
and 95% (light) confidence intervals of our MCMC fit. The fact that such an agreement can be achieved is an important success given 
the simplicity of the formation history model. 




logio(Mpeak/M0) 



# of propositions 



> 1388 




logio(Kiax/kms \ 



# of propositions 



1050 1200 1350 



1050 1200 1350 



Figure 6. Marginalised posterior probability distributions for the Mvir (left) and Vmax (right) static model parameters when calibrating 
against the red and blue stellar mass functions at z=Q (Fig. 5). All panels have been zoomed in to the high probability regions. Contours 
on the 2-D (blue) panels indicate the 68 and 95% confidence regions. Yellow dots mark the marginalised most probable parameter values. 
The diagonal panels show the marginalised 1-D distributions with the 68 and 95% confidence intervals (dark and light shaded regions 
respectively). The approximately Gaussian shape of these 1-D distributions indicates the well behaved nature of the model. The only 
parameter degeneracies are between the normalisation (^M^ir /^Vmax ) ^^^ width (c^M^ir /<^Vmax ) of the physics functions. 



at approximately the correct position in colour space (Cole 
et al. 2005) is an interesting result for such a simple model. 
In the left hand panel of Fig. 5 we show the red and blue 
model galaxy stellar mass functions (solid lines) against the 
corresponding constraining observations (error-bars). De- 
spite its simplicity, the chosen form of the physics function 
produces a good reproduction of the data. This is true across 
a wide range in stellar mass, indicating that the model is 
capable of successfully matching the integrated time evolu- 
tion of stellar mass growth as a function of halo mass at 
2;=0. Also, since blue galaxies preferentially trace those ob- 
jects which have undergone significant recent star formation, 
the model's reproduction of the observed blue mass func- 



tion suggests that the rate of star formation as a function of 
stellar mass near z—0 is also in broad agreement with the 
observed universe. 

The fact that such an agreement is attainable with 
this simple model should be viewed as a key success of the 
methodology and a validation of the general form we have 
chosen for the physics function, Fphys(Mvir). Having said 
this, there are some differences in the left hand panel of 
Fig. 5 worth noting. In particular, there is an over-prediction 
in the number density of the most massive red galaxies and 
a corresponding under prediction of the most massive blue 
galaxies. 
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3A.2 The Vm 



model 



Having established that a physics function constructed using 
Mvir as the single input variable can successfully provide a 
good match to the observed z=0 red and blue stellar mass 
functions, we now turn our attention to the results of using 
Vmax as the input property (Eqn. 4). 

In the right hand panel of Fig. 5 we present the colour 
split stellar mass functions for the static Knax model. Again, 
a fixed colour division of B—V=0.8 is used to define the two 
colour populations and we use MCMC tools to constrain 
the physics function parameter values to provide the best 
statistical reproduction of the Bell et al. (2003) data. The 
resulting parameter values are presented in Table 1. Un- 
surprisingly, a comparison with the equivalent values of the 
Mvir model indicates that the peak efficiency of convert- 
ing fresh baryonic material into stars in a single time step 
remains similar (^Vm^^=0.53; Sm^^^^O.^G). However, the av- 
erage virial mass of halos with ymax~(lpeak=158 kms~^) is 
10^^'^ M0, therefore this peak efficiency occurs in slightly 
more massive halos than was the case for the Mvir model. 
This is a reflection of the different growth histories of these 
two halo properties. 

As was the case for the Mvir model, an excellent repro- 
duction of the observations is attainable when using Vmax as 
the input parameter to the physics function. We find that 
the over-prediction of high mass red galaxies has been allevi- 
ated, although at the cost of now somewhat under-predicting 
the number density of low mass blue galaxies. Importantly 
though, given a suitable choice for values of the free parame- 
ters of the physics function, both the Mvir and Vmax physics 
functions can produce a good match to the distribution and 
late time growth of stellar mass at z=0 despite the differ- 
ences in their mean time evolution (cf. Fig. 3). 

In Fig. 6 we present the marginalised posterior proba- 
bility distributions for our MCMC calibration of both the 
Mvir (left panel) and Vmax (right panel) models. For clar- 
ity we have zoomed in on the regions of high probability 
in all panels instead of showing the full ranges explored. 
The well behaved and understandable nature of the param- 
eter distributions gives us further confidence in the validity 
of our model implementation. The approximately Gaussian 
shape of the 1-D distributions (diagonal panels), coupled 
with their uni-modal nature, indicates that all of the pa- 
rameters are well constrained. Furthermore, the 2-D pan- 
els demonstrate that the only degeneracies in either model 
are between those parameters controlling the normalisation 

(^Mvir /^Vmax ) ^^^ wldth (cTM^ir /^Vmax ) ^^ ^^C log-UOrmal 

physics function. This makes intuitive sense as these param- 
eters jointly determine the integral of the star formation 
rate defined by Eqn. 1 and therefore the approximate total 
amount of stellar mass formed by each galaxy. 



3.2 High redshift 

In the previous section we demonstrated that our simple 
formation history model is capable of reproducing the ob- 
served red and blue galaxy stellar mass functions of the lo- 
cal Universe. We also showed that this is true independent 
of whether we utilise the Mvir or Knax form of the physics 
function (Eqns. 3 & 4). However, as seen in Fig. 3, there 
are important differences in the time evolution of these halo 
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Figure 7. The 2;~0.9, 1.8, 3.3 stellar mass functions predicted 
by the static Mvir (dashed black line) and Vmax (dashed orange 
line) models (Eqns. 3,4). Observational data from Perez-Gonzalez 
et al. (2008, PG08) and is shown for comparison. The solid lines 
indicate the results of convolving the model stellar masses with a 
normally distributed random uncertainty of 0.3 or 0.45 dex (for 
redshifts less than/greater than 3 respectively) in order to mimic 
the systematic uncertainties associated with the observed masses. 
The dark and light shaded regions show the 68 and 95% confi- 
dence intervals predicted using the marginalised parameter dis- 
tributions obtained from our z=0 parameter calibration (Fig. 6). 
There are clear differences between the Mvir and Vmax models 
at higher redshifts, despite their approximate agreement at z=0 
(Fig. 5). This reflects the different time evolution of these two 
halo properties (see Fig. 3). 



properties. This suggests that we should see corresponding 
differences in the galaxy populations predicted at higher red- 
shifts. 

In Fig. 7 we present the stellar mass functions of both 
the Mvir and Vmax models (dashed lines) against the ob- 
served z>0 mass functions of Perez-Gonzalez et al. (2008, 
points). The solid lines represent the formation history 
model results after a convolution with a normally distributed 
random error of dispersion 0.3 dex for z<3 and 0.45 dex for 
redshifts greater than this value (Moster et al. 2013). Such a 
convolution is common practice and approximates the miss- 
ing uncertainties in the observational data due to systemat- 
ics involved with producing stellar mass estimates from high 
redshift galaxy observables (e.g. Fontanot et al. 2009; Guo 
et al. 2011; Santini et al. 2012). 

There are clear quantitative differences between the 
stellar mass functions produced by the two models. These 
become more pronounced as we move to higher redshifts. At 
z^3 (bottom panel), the Vmax model predicts a sharp fall off 
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Figure 8. The mean evolution of the stehar-halo mass relation of central galaxies. The results of the static Mvir (Eqn. 3; black solid line) 
and Vmax (Eqn. 4; orange dashed line) models are both shown. The dark and light shaded regions indicate the 68 and 95% confidence 
intervals predicted using the marginalised parameter distributions from our calibration at z=^ (see Fig. 6). A comparison with the 
subhalo abundance matching results of Moster et al. (2013) (blue error-bars) indicates that both forms of the physics function fail to 
reproduce the evolution of the stellar mass growth efficiency required to reproduce the high redshift stellar mass function (see Fig. 7). 



in the number density of galaxies with stellar masses greater 
than 10^°"^ Mq. When using Mvir to define the physics func- 
tion, this drop off does not occur until M^^^IO^^M©, re- 
sulting in a differing prediction in the number density of 
these galaxies by greater than 2 orders of magnitude at high 
masses. Despite this, both versions of the physics function 
predict 2;>0 stellar mass functions which are too steep at 
high masses, although the addition of the random uncer- 
tainties (solid lines) largely alleviates this problem. There 
are also notable differences at lower stellar masses, where 
both models over-predict the number of galaxies. The Vmax 
model also predicts that a large fraction of galaxies with 
Mvir < 10^°"^ Mq are already in place by z—?>^ with a corre- 
spondingly slower evolution to z—0. 

Many of these differing qualitative predictions can be 
understood by considering the differences in the time evo- 
lution of Mvir and Vmax as shown in Fig. 3. For example, 
the deficit of high stellar mass galaxies in the Vmax model at 
z^?t is due to their host halos initially being identified with 
Vmax values greater than Vpeak (at least for the mass reso- 
lution of our simulation) . These halos therefore spend little 
time in the efficient star forming band. The result is a re- 
duced amount of in-situ star formation at early times, with 
effects that carry all the way through to ^=0 as these galax- 
ies grow, predominantly through merging. We can similarly 
understand the cause of the larger predicted number density 
of high redshift low mass galaxies in the Vmax model. In this 
case, the lowest mass halos present at high redshift s have 
spent a longer time close to the peak of the efficient star 
forming band. This results in these halos already hosting 
significant amounts of stellar mass by z—?>. 

To illustrate this further, in Fig. 8 we show the evolution 
of the mean stellar-halo mass relation for both models. The 
blue error-bars represent the relations predicted by the sub- 
halo abundance matching model of Moster et al. (2013). We 
have specifically chosen to compare our results against the 
work of Moster et al. (2013), as they take their halo masses 
from the same dark matter merger trees as used in this work 
(as well as the higher resolution Millennium-II simulation; 
Boylan-Kolchin et al. 2009) and also construct their model 
to match the same high-redshift stellar mass functions of 
Perez- Gonzalez et al. (2008). Hence, the blue error-bars of 
Fig. 8 represent the evolution in the integrated stellar mass 
growth efficiency which our model must achieve in order to 



successfully replicate the observed stellar mass functions of 
Fig. 7. 

By construction, both the Mvir and Vmax models pro- 
duce extremely similar relations at z—0 but with clear differ- 
ences at higher redshifts. It is these variations in the typical 
amount of stars formed within halos of a given mass that 
drives the different predictions for the evolution of the stel- 
lar mass function. For example, the much higher average 
stellar mass content of low mass halos at z=3 when using 
the Vmax model (Fig. 8) is the cause of the increased nor- 
malisation of the low mass end of the relevant stellar mass 
function in Fig. 7. 

Importantly, it can be seen from Fig. 8 that neither 
the Mvir nor the Vmax model reproduce the evolution of the 
stellar-halo mass relation found by Moster et al. (2013); in 
particular the position and normalisation of the peak value. 
The use of a redshift independent halo mass to define the 
peak in-situ star formation efficiency of the Mvir model re- 
sults in no change to the position of the peak of the stellar- 
halo mass relation with redshift. Although not immediately 
obvious why this should be so, it can be understood by con- 
sidering the typical evolution of a halo across the relation. 

At early times, halos grow in mass rapidly, however 
they typically still sit below the efficient star formation mass 
regime defined by the physics function (see Fig. 3). In Fig. 8, 
these halos will therefore travel almost horizontally from left 
to right with a low stellar-halo mass fraction. Eventually 
halos will enter the mass regime of efficient star formation, 
causing them to rapidly increase their stellar-halo mass frac- 
tions with only a relatively modest growth in halo mass. 
This phase of rapid stellar mass growth causes a 'pile-up' of 
galaxies in the stellar-halo mass relation that peaks around 
the virial mass at which halos again transition out of the 
efficient star forming regime. Since the mass at which this 
occurs is fixed in our simple static model, the position of 
the stellar-halo mass relation peak is therefore also fixed in 
the Mvir model. Due to the evolving Mvir-Vmax relationship, 
the position of the peak efficiency for the Vmax model does 
evolve, but unfortunately in the direction opposite to that 
required. The shallower tail of the relation towards higher 
halo masses is caused by the subsequent growth of galaxies 
due to mergers. 

Based purely on this inability to reproduce the required 
evolution in the stellar-halo mass relation, it is clear that 
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Figure 9. Marginalised posterior probability distributions for the Mvir (left) and Vmax (right) parameters of the redshift dependent 
physics function. In both panels, the models were constrained to simultaneously reproduce the observed z=0 red and blue stellar mass 
functions (Fig. 12) as well as the time evolution of the stellar-halo mass relation (Fig. 9). All panels have been zoomed in to the high 
probability regions. Contours on the 2-D (blue) panels indicate the 68 and 95% confidence regions. Yellow dots mark the marginalised 
most probable parameter values. The diagonal panels show the marginalised 1-D distributions with the 1 and 2-a confidence intervals 
shown by light and dark shaded regions respectively. The approximately Gaussian shape of these 1-D distributions indicates the well 
behaved nature of the model. 



there is no chance for our current, non-evolving physics func- 
tion to match the observed distribution of stellar masses in 
both the low and high redshift Universe simultaneously, irre- 
spective of the values of the available parameters or whether 
Mvir or Vmax IS uscd as the dependant variable. 

3.3 Incorporating a redshift evolution 

Although capable of reproducing the observed red and blue 
stellar mass functions at ^=0, we showed in §3.2 that our 
simple formation history model struggles to reproduce the 
high redshift distribution of stellar masses. Importantly, we 
also concluded that there is no combination of physics func- 
tion parameter values (see Eqns. 3 & 4) which can alleviate 
this discrepancy. In this section we therefore look to extend 
our simple model by introducing a redshift dependence to 
the physics function. In effect, this is equivalent to the intro- 
duction of an evolution of the star formation efficiency with 
time for a fixed halo mass/maximum circular velocity. Such 
an evolution is well motivated both theoretically and obser- 
vationally, suggesting the presence of alternative/additional 
star formation mechanisms at high redshift when compared 
to those of the local Universe. For example, so called "cold- 
mode" accretion (Birnboim & Dekel 2003; Keres et al. 2005; 
Brooks et al. 2009) is thought to be able to effectively fuel 
galaxies of massive halos at high redshift, allowing for in- 
creased star formation. In addition, the early universe was 
a more dynamic place, with an enhanced prevalence of gas 
rich galaxy mergers and turbulence driven star formation 
(e.g. Dekel et al. 2009; Wisnioski et al. 2011). 

To reproduce the evolving position and normalisation 
of the stellar-halo mass relation as found by Moster et al. 
(2013), we modify the physics function of Eqn. 3 by intro- 
ducing a simple power law dependence on redshift to each 
of the free parameters: 

(6) 






(7) 
(8) 



logio(Mpeak(2;)) = logio(Mpeak) " (1 + 2^) 



«M,, 



where at z=0: Mpeak(2:)=Mpeak, cFM^;^{z)=aM^;^ and 
^Mvir('2^)=^Mvir- The exact values of the redshift scalings 
are calibrated using MCMC to provide the best simulta- 
neous reproduction of the Moster et al. (2013) stellar-halo 
mass relation at 2;=0, 1, 2 & 3, as well as the ^=0 red and 
blue stellar mass functions, and are presented in Table 1. 

In the left hand panel of Fig. 9 we present the relevant 
marginalised posterior probability distributions of the six 
free model parameters. Similarly to the redshift-independent 
case (cf. §3.1), the approximately Gaussian shape of the 1-D 
probability distributions indicates that the parameters are 
generally well constrained. However, in addition to the de- 
generacy between the physics function normalisation (^M^ir ) 
and width (aM^ir) noted in §3.1, there are also clear and 
understandable degeneracies between the redshift evolution 
and z—0 value of each parameter (e.g. Mpeak and aM^ir)- 

We note that there are minor differences between the 
z—0 stellar mass function utilised by Moster et al. (2013) to 
constrain their stellar-halo mass relation (Li & White 2009), 
and the z—0 mass function which we employ in this work 
(Bell et al. 2003). However, we calibrate our model against 
both the stellar-halo mass relation and colour-split stel- 
lar mass functions at z—0 with equal weights. The MCMC 
fitting procedure then attempts to find the best compro- 
mise between these two (as well as all other) constraints. 
Since we find that there are no multi-modal features in the 
marginalised posterior probability distributions (see Fig. 9), 
the parameter sets required to fit each constraint individ- 
ually must be statistically compatible with each other. We 
therefore conclude that this slight inconsistency in our fit- 
ting procedure has minimal effect on our ability to demon- 
strate the success and utility of the model and on our results. 

The preferred values of o^M^ir and /^M^ir ^^^ relatively 
small, indicating little need for evolution in both the peak 
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Figure 10. The evolution of the stellar-halo mass relation of central galaxies for the evolving Mvir (Eqn. 3; black solid line) and T/max 
(Eqn. 4; grey dashed line) models. Orange shaded regions indicated the subhalo abundance matching results of Moster et al. (2013). 
A comparison with Fig. 8 indicates that by shifting both the normalisation (^M^ir^ ^Vmax) ^^^ position (Mpeak? ^peak) of the physics 
functions, we are able to reproduce the correct evolution of the stellar-halo mass relation at high redshifts in terms of both shape and 
amplitude. This leads to a much better agreement between the observed and predicted stellar mass functions at z>Q. (see Fig. 11). 



position, Mpeak(2;), and width, (tm^-^^{z)^ of the physics func- 
tion. As a consequence, the values of both Mpeak and cTM^ir 
are similar to the non-evolving case (cf. Table 1). How- 
ever, there is a strong evolution preferred for the normal- 
isation of the physics function, ^m^^^^ such that it decreases 
rapidly with increasing redshift. In order to maintain the 
total 2;=0 stellar mass density, the value of Mpeak=0.9 is 
therefore considerably higher than was the case in the non- 
evolving form of the physics function. This implies that 
90% of all freshly accreted baryonic material in halos with 
log ]^Q (Mpeak 7^0) = ! 1.4 is converted into stars at z—0. How- 
ever, at z—1^2 and 3 the peak conversion efficiencies are con- 
siderably lower: 54%, 40% and 32% respectively. 

We also similarly modify the Vmax physics function, 

-^phys V "^ max j • 



loglo(^peak(2)) 



Iogio(ypeak)-(1+;Z)"^— , (9) 

<TV„ax • (1 + ZY""--^ , (10) 

^V„.. • (1 + ^r"— , (11) 



with the redshift scalings being calibrated to reproduce the 
same observations as the Mvir case above. The marginalised 
posterior probability distributions are presented in the right 
hand panel of Fig. 9, with the preferred parameter values 
again presented in Table 1. 

As was found for the redshift dependent Mvir model, 
the marginalised posterior distributions indicate that the 
free model parameters are well constrained and that there 
are no unexpected degeneracies between them (Fig. 9). If we 
consider the preferred values of the parameters themselves, 
we again find that the normalisation of the physics function, 
^Vmax(^)7 shows the most pronounced evolution. A value of 
7Vmax— ~ 0-98 indicates that the maximum star formation 
efficiency declines almost linearly as a function of (l-\-z). 
This strong evolution requires a value of ^v^ax ^^ ^—^ that 
is actually greater than 1, and hence more than the total 
freshly accreted baryonic material must be converted into 
stars in halos with log;Lo(^peak/(kms~"^))=2.1 at this red- 
shift. This suggests that we may need to vary the universal 
baryon fraction as a function of halo mass (or maximum 
circular velocity), perhaps to mimic the effects of processes 
such as the recycling of ejected baryons during star forma- 
tion (e.g. Papastergis et al. 2012). 

In Fig. 10, we present the stellar-halo mass relations 
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Figure 11. The 2;~0.9, 1.8, and 3.3 stellar mass functions pre- 
dicted by the evolving Mvir (dashed black line) and Vmax (dashed 
orange line) models. Observational data from Perez-Gonzalez 
et al. (2008) is shown for comparison. The solid lines give the 
results of convolving the model stellar masses with a normally 
distributed random uncertainty of 0.3 or 0.45 dex (for redshifts 
less than/greater than 0.3 respectively) in order to mimic the 
systematic uncertainties associated with the observed masses. By 
using an appropriate redshift evolution of the physics function pa- 
rameters, the model's ability to successfully recover the observed 
high redshift stellar mass functions is improved. 



of the new, redshift dependant, Mvir (black) and Vmax (or- 
ange) models. The blue error-bars again indicate the results 
of Moster et al. (2013). By incorporating the redshift de- 
pendence we are now able to successfully reproduce the evo- 
lution of both the normalisation and peak position of the 
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Figure 12. The red and blue galaxy stellar mass functions produced by the evolving Mvir (left) and T/max (right) formation history 
models. Error-bars indicate the observations of Bell et al. (2003). The free model parameters have been constrained to simultaneously 
reproduce these mass functions as well as the evolution of the stellar-halo mass relation (Fig. 10). The dark and light shaded regions 
show the associated 68 and 95% confidence regions obtained from this calibration. 



stellar-halo mass relation required at z^O. The effects of 
this on the predicted high redshift stellar mass functions 
of both the Mvir and Vmax models can be seen in Fig. 11. 
As expected, we now find an improved agreement with the 
observations when compared to the original, non-evolving 
physics function results (cf. Fig. 7). 

For completeness we also present the z=0 colour-split 
stellar mass functions for both models in Fig. 12. A reason- 
able agreement with the constraining observational data is 
still achieved. However, we now find that an under prediction 
in the number density of massive blue galaxies is present in 
both models, suggesting that our implemented evolutionary 
model may not provide enough late time star formation in 
the most massive halos. 



4 DISCUSSION 

In §3.1, we demonstrated that our most basic, non-evolving 
form of the physics function is able to successfully reproduce 
the observed red and blue stellar mass functions of the local 
Universe (Fig. 5). This key result highlights the utility and 
validity of our basic methodology and model implementa- 
tion. In addition, it reinforces the commonly held belief that 
the growth of galaxies is intrinsically linked to the growth 
of their host dark matter halos (White & Rees 1978). 

Although the level of agreement achieved with the ob- 
served z=0 colour split stellar mass functions is generally 
very good, there are some discrepancies. For example, there 
is an under-prediction in the number density of the most 
massive blue galaxies in the Mvir model (left panel of Fig. 5), 
with a corresponding over-prediction in the number of mas- 
sive red galaxies. Our z>0 analysis suggests that this is at 
least partially due to an incorrect evolution of the stellar- 
halo mass relation with time (see Fig. 8). However, we also 
note that an excess of massive red galaxies is a common fea- 
ture of traditional semi-analytic galaxy formation models 
which similarly tie the evolution of galaxies to the masses of 
their host dark matter halos. In such models, efficient feed- 



back from active galactic nuclei is typically responsible for 
truncating star formation in the most massive galaxies and 
hence causes the average stellar populations of these objects 
to become older and redder (Croton et al. 2006; Bower et al. 
2006; Mutch et al. 2011). This is already mimicked within 
the framework of our formation history model through the 
turnover at the high Mvir (or Vmax) end of the physics func- 
tion. However, a more gradual cutoff may be required in the 
Mvir model case, in order to allow star formation to proceed 
for longer in the galaxies populating the most massive halos. 

In this paper we have deliberately restricted ourselves 
to considering only a very simple form of the physics func- 
tion. This has allowed us to take advantage of the resulting 
transparency when interpreting our findings. However, we 
stress that the model can easily be extended to include ar- 
bitrary levels of complexity. For example, we have chosen 
to use a log-normal distribution to define the form of the 
physics function. Although being conceptually simple, the 
symmetric nature of this formalism implicitly assumes that 
the physical mechanisms responsible for quenching star for- 
mation in both low and high mass halos scale identically 
with halo mass (Fig. 3). This assumption has little physical 
justification and in order to provide the best results, it may 
be necessary to independently adjust the slope of the func- 
tion at both low and high masses, and perhaps even as a 
function of redshift. In future work, we will address this is- 
sue by carrying out a full statistical analysis aimed at testing 
a number of different functional forms for both the physics 
and baryonic growth functions. 

However, even within the reduced scope of this current 
work, we have learned a great deal from simply examining 
the high redshift stellar mass function predictions of the 
formation history model. In particular, we have highlighted 
the need for the physics function to produce an evolution 
in the stellar-halo mass relation as a function of redshift 
in order to match the observed space density of massive 
galaxies at early times. Using Vmax as the input property to 
the function introduces such an evolution, but in the wrong 
direction. Future improvements to the model could focus on 
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finding a iialo property tiiat does evolve correctly with time 
and would thus be a more natural anchor of the physics 
function. This would avoid the need to artificially introduce 
an evolution to match the observations, as we have done 
here. 

Although the need for an evolving stellar-halo mass re- 
lation has been discussed in the literature, the precise form 
with which this evolution manifests itself is less clear. The 
results of subhalo abundance matching studies, such as that 
of Moster et al. (2013) (which we compare to in this work), 
are quite sensitive to the choice of input datasets and the 
technical aspects of the methodology. For example, Moster 
et al. (2013) find that the peak stellar-halo mass ratio in- 
creases from just 0.15% at z =4 to 4% at ^=0, with a corre- 
sponding shift in position from a halo mass of 10^^"^ M© to 
10^^'^ M0. In contrast, an alternative study carried out by 
Behroozi et al. (2012) finds very little change in either the 
normalisation or peak location over a broad range in red- 
shift. However, they do note a marked drop in the relation 
for the most massive halos at 2;<2.5. This results in a qualita- 
tively different prediction for the evolution of these massive 
halos, such that their efficiency of converting baryons in to 
stars is higher as a function of look-back time (the opposite 
trend to that found by Moster et al. (2013)). 

A potentially valuable use for the formation history 
model is to provide a general consistency check of subhalo 
abundance matching studies. The physics function could 
be adapted to exactly replicate the shape and evolution 
of the star formation efficiencies they predict (Behroozi 
et al. 2013), allowing their validity to be assessed when self- 
consistent ly applied to individual dark matter merger trees. 
The additional galaxy properties provided by the formation 
history model, such as star formation histories and colours, 
could be used to further compare and contrast the success- 
fulness of different abundance matching methodologies. For 
example, it has been suggested that galaxy mergers may re- 
sult in a significant fraction of the in-falling satellite stellar 
mass being added to a diffuse intra-cluster light component 
instead of to the newly formed galaxy (Monaco et al. 2006; 
Conroy et al. 2007). The strength of this effect is expected 
to increase significantly with increasing halo mass and is in- 
cluded in the subhalo abundance matching study of Behroozi 
et al. (2012) but not Moster et al. (2013). By simply adding 
a mass dependent amount of stellar material to an ICL com- 
ponent during merger events, our simple formation history 
model could be easily adapted to explore such a scenario. 

We also note that the simplicity of our formation history 
model results in it being extremely fast and computation- 
ally inexpensive when compared to traditional semi-analytic 
models. This has allowed us to straightforwardly calibrate it 
against a number of observed relations using Markov Chain 
Monte-Carlo techniques. This procedure can also be trivially 
extended to provide statistically accurate (against select ob- 
servations) mock catalogues for use with large surveys. Fur- 
ther to what can be achieved using current Halo Occupation 
Distribution or sub-halo abundance matching methods, cat- 
alogues produced using our model include both full growth 
histories and star formation rate information for each indi- 
vidual galaxy, with no need to add any artificial scatter to 
approximate variations in formation histories. In addition, 
the direct and clear dependence of the model on the halo 
properties of the input dark matter merger trees makes it 



an ideal tool for investigating a number of additional top- 
ics. Examples include: comparing the effects of variations 
between different N-body simulations and halo finders on 
the physics of galaxy formation and evolution; investigating 
the predictions of simple monolithic collapse scenarios; con- 
trasting various mass dependent merger star burst models; 
and exploring the ramifications of N-body simulations run 
with alternative theories of gravity. 

Finally, we note that a common criticism of semi- 
analytic models is the presence complex degeneracies be- 
tween large numbers of free parameters. The MCMC cal- 
ibration procedure we employ highlights the complete ab- 
sence of such degeneracies in our formation history model, 
further demonstrating its well behaved and understandable 
nature. 



5 CONCLUSIONS 

In this work we introduce a simple model for self-consistently 
connecting the growth of galaxies to the formation history 
of their host dark matter halos. This is achieved by directly 
tying the time averaged change in mass of a halo to the 
star formation rate of its galaxy via two simple functions: 
the "baryonic growth function" and the "physics function" 
(Eqns. 2,3). We utilise N-body dark matter merger trees 
to provide self consistent growth histories of individual ha- 
los that naturally includes scatter due to varying formation 
histories. This allows us to produce full star formation histo- 
ries for individual objects, and thus provide predictions for 
secondary properties such as galaxy colour. 

While closely related to other models in terms of its ba- 
sic methodology (Bouche et al. 2010; Cattaneo et al. 2011), 
our model has a number of important generalisations which 
enhance its utility. In particular, we implement a single, uni- 
fied physics function which encapsulates the effects of all of 
the intertwined baryonic processes associated with galactic 
star formation and condenses them down into a simple map- 
ping between star formation efficiency and dark matter halo 
properties. The qualitative form of this function is motivated 
by our general knowledge of galaxy evolution, however, in 
this work we make no attempt to directly tie it to individ- 
ual physical processes or their particular scalings with halo 
properties. 

As well as introducing this new model, we demonstrate 
its ability to replicate important observed relations such as 
the galactic stellar mass function, and also illustrate some 
examples of its potential for investigating different theories 
of galaxy formation and evolution. Our main results can be 
summarised as follows: 

• Motivated by the observed suppression of star forma- 
tion efficiency in both the most massive and least massive 
dark matter halos we begin by parametrising the physics 
function as a simple, non-evolving, log-normal distribution 
with a single independent variable of either halo virial mass, 
Mvir, or maximum circular velocity, Vmax (Fig. 3). 

• With just three free parameters controlling the posi- 
tion, normalisation and dispersion of the peak star forma- 
tion efficiency, we show that the formation history model 
can successfully reproduce the observed red and blue stellar 
mass functions at redshift zero. Assuming a suitable choice 
of the parameters, this result is independent of the use of 
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Mvir or Vmax as the dependant variable of the physics func- 
tion (Fig. 5). 

• For the purposes of repUcating the stellar mass func- 
tions across a wide range of redshifts, we find our static 
model to be inadequate. This is due to its inability to pro- 
duce the correct evolution of the stellar-halo mass relation 
with time (Figs. 7 & 8). 

• We therefore investigate the use of redshift as a sec- 
ond dependant variable to the physics function in order to 
control the position and normalisation of the peak star for- 
mation efficiency with time. Using this simple adaptation 
alone, the formation history model is able to better repro- 
duce the observed high redshift stellar mass functions out 
to ^=3.5 (Figs. 10 & 11) whilst still maintaining a good 
reproduction of the z=0 colour-split stellar mass function. 

• By statistically calibrating the free model parameters 
using Markov Chain Monte-Carlo chain techniques through- 
out this work, we are able to use the marginalised posterior 
likelihood distributions to demonstrate the well behaved and 
transparent nature of our simple model (Fig. 6). 

In order to demonstrate its construction and utility we 
have presented one of the simplest forms of the formation 
history model. However, a fundamental strength of its con- 
struction is that it can be easily extended to arbitrary levels 
of complexity in order to investigate a whole host of physical 
processes associated with galaxy formation and evolution, 
some general examples of which we have outlined in §4. In 
future work we will investigate the predictions made when 
using alternative forms of the baryonic growth and physics 
functions. We will also extend the model to investigate the 
birth of super-massive black holes and the evolution of the 
quasar luminosity function. 
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